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Given an unconditionally stable algorithm for solving the Cahn-Hilliard equation, we present 
a general calculation for an analytic time step Ar in terms of an algorithmic time step At. By 
studying the accumulative multi-step error in Fourier space and controlling the error with arbitrary 
accuracy, we determine an improved driving scheme At — At^^^ and confirm the numerical results 
observed in a previous study 
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The Cahn-Hilliard (CH) equation H models the phase 
separation that occurs during the quench of a conserved 
system from a high temperature isotropic phase into two 
distinct phases at low temperatures. The pattern of 
the two phase regions coarsens as the time r increases, 
i.e. the length-scale of these regions grows. At the 
later stages of this phase ordering process, the dynam- 
ics are dominated by a single length scale, the pattern 
domain size L, which increases with a power law in time 
T, L{t) ^ T^/^ 2\. This power-law growth implies that 
the motion of the domain walls becomes extremely slow 
at late times after a quench since a typical domain wall 
speed is V ^ dL / dr ^ r^^/"^, and a typical time scale for 
the interface to move a distance of order the interfacial 
width ^ is of order ^/v ^ t^I'^ . 

Since there is no known analytic solution of the Cahn- 
Hilliard equation for random initial conditions, computa- 
tional methods are necessary for investigation of such sys- 
tems. The most straightforward approach is the Euler al- 
gorithm, which must employ a time step ^Ieu ^ (Ax)'* if 
stability is to be maintained, where Aa; is the lattice spac- 
ing. Additionally, for Cahn-Hilliard systems, to resolve 
the interfacial profile, one has to use a lattice spacing 
Aa; < The Euler fixed time step is suitable for update 
near the interface but wastefuUy accurate in the bulk at 
late times. This has been the main challenge of computer 
simulation of Cahn-Hilliard systems. The recently devel- 
oped unconditionally stable algorithm Q |^ elegantly 
overcomes this difficulty. It allows a mode- dependent ef- 
fective time step Atg/ / — a larger effective time step in 
the bulk as the domain size gets larger, while keep the 
effective time step finite near the interfacial region (see 
Eq. (Q). Since the unconditionally stable algorithm al- 
lows no constraints on time step, the main issue is ensur- 
ing the accuracy of the simulation. In a previous study 
[5j , Cheng and Rutenberg numerically demonstrated that 
the error in correlations decreases monotonically as A de- 
creases down to A = 0.001, where ^ is a prefactor in 
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the previous driving scheme At — Atg and tg is the 
structural time (see below for a precise definition) . How- 
ever, absence of computational power prevents us from 
exploring the error behavior for arbitrarily small A. For 
arbitrary accuracy, we need to rigorously prove that the 



error can be made arbitrarily small. 

To achieve this goal, we must first distinguish two 
quantities generic to all numerical algorithms: analytic 
time step Ar [analytic time r) and algorithmic time step 
Ai [algorithmic time t). The former appears in the equa- 
tion of motion and represents the time step (time) of the 
system evolution governed by the exact solution to the 
dynamical equations, while the latter appears in the finite 
difference scheme and represents the time step (time) of 
the system evolution by the computational algorithms. 
We now briefiy review these two concepts and study why 
the distinction has been largely overlooked thus far. 

In Euler algorithm, it is not necessary to distinguish 
the analytic time step and the algorithmic time step, 
since there is a threshold on the time step, and they 
are always approximately identical (see the analysis be- 
low). On the other hand, semi- implicit algorithms have 
no such threshold. An unconditionally stable algorithm, 
an extension of the semi-implicit method, allows for an 
arbitrarily large time step without encountering the nu- 
merical instabilities for suitably chosen parameters (as 
determined using a standard von Neumann stability anal- 
ysis). Although this method has been in use for some 
time, there has been little analytic study about how to 
obtain maximal speedup while controlling the accuracy. 
Indeed, all the previous work known to us assumes no 
difference between the analytic time step and the algo- 
rithmic time step — one can only increase the time step 
modestly and assume the resulting error are small enough 
to be ignored. 

In what follows, we perform a general calculation of the 
analytic time step Ar in terms of the algorithmic time 
step At, and show how this relationship allows one to 
choose a driving scheme for arbitrary accuracy. Concomi- 
tantly, we demonstrate that the driving scheme can be 
improved to At — At^^^ . While the calculation presented 
is specifically applicable to the Cahn-Hilliard equation, 
much of our analysis is general, and should guide subse- 
quent studies of more complicated systems. For simplic- 
ity but without loss of generality, we restrict our analysis 
to two dimensions (2D). 
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The Cahn-Hilliard equation can be written as 



where the free energy functional is 



F= i (fx 
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If 
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and 0(x, r) is a conserved scalar field (such as an appro- 
priately scaled mass concentration) and the potential has 
a double- well structure that the equilibrium values are at 
(/) = ±1. To illustrate and distinguish the analytic time 
step from the algorithmic time step, we first study the 
exact dynamics of the Cahn-Hilliard systems Eq. |^ in 
Fourier space. At an analytic time t, the system evolu- 
tion after an analytic time step Ar is governed by the 
Taylor expansion: 
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Using the results of the field derivatives in Fourier space 
Eq. (|15|l . for a finite Ar, one finds that all n > 2 terms 
are negligible compared with the n = 1 term. So we 
obtain the traditional Eulcr finite difference scheme 



(j)k{t + AtEu) = Mt) + ^^Eu^ 

OT 



(4) 



where d(f)k/dT is the Fourier transform of d(t>/dT in Eq. 

and is a function of 4>k{t)- We see that the Euler al- 
gorithm uses first order finite differences to approximate 
the solution obtained by exact dynamics. On the other 
hand, an unconditionally stable algorithm is obtained by 
an appropriate semi-implicit discretization of Eq. in 
algorithmic time: 

^t+At + (1 - ai)AtV'^(j)t+At + (1 - a2)AtV^(t)t+At 

^(j)t- AtV^{ai,Pt+a2V^<pt~<ft)- (5) 

Unconditionally stability is obtained for the choices ai > 
2 and 02 < 0.5 0. Here, (t)t+/\t represents the implicit 
terms and 4>t represents the explicit terms. We can solve 
Eq. (O directly in Fourier space and obtain 



0fc(i + At) = c^k{t) + At,ff{k, At)- 



dT ' 



where the fc-dependent effective time step is 



Ate//(fc,At) 



At 



l-AtAk[(ai-l) + (a2-l)Ak]' 
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and Ak = —k^ is the Fourier-transformed Laplacian. The 
Euler algorithm has a mode-independent fixed time step 
to update the system in Fourier space, but, as Eq. 
reveals, the unconditionally stable algorithm has a mode- 
dependent effective time step Ateff{k,At). A direct 



comparison of Eq. (0 and Eq. yields that the ana- 
lytic time step Ar is always a good approximation of the 
algorithmic time step Atsu in Euler algorithm. However, 
for the unconditionally stable algorithm, a comparison of 
Eq. ^ and Eq. ^ docs not give a straightforward rela- 
tion between Ar and At, i.e., we do not know what Ar 
corresponds to At. In what follows we explore the rela- 
tionship between these two time steps in CH equation, 
and the consequences this relationship has on the accu- 
racy of the solution method. The steps of our procedure 
shown in italics. 

Calculate the analytic time and time step: We now 
calculate the analytic time step Ar in terms of an algo- 
rithmic time step At. Cahn-Hilliard systems are purely 
dissipative systems — the energy density E monotoni- 
cally decreases with the analytic time with the relation 
E cx r~^/^ 2]. Without such a relationship between a 
physical quantity and the analytic time, the analysis that 
is performed below cannot proceed, and thus progress in 
applying these methods to other models hinges on the 
physical insights needed to obtain such relationships (in 
this case the so-called "scaling hypothesis"). The ana- 
lytic time is conveniently calculated in terms of the mono- 
tonically decaying energy density E: t = B/E^, where 
the prefactor B can be numerically determined by requir- 
ing Ar = At as At — !■ in the late-time scaling regime 
since our unconditionally stable algorithm is arbitrarily 
accurate as At 0. Note that the calculation here is 
identical to the calculation of the structural time tg in a 
previous study |5j since the structural time is just another 
representation of the analytic time. 

We can calculate the analytic time step by differenti- 
ating r with respect to E: 



AE r4/3 
^---^B^^-3AE^^, 
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and AE can be calculated by integrating AE from each 
Fourier mode: 



AE 
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;Ateff{k,At)Tk, (9) 



where the time derivative dcji-k/dr = —k'^SF/S<pk from 
Eq. (P) and Acjjk = Mt + At) - (j^tit) = Ateffdcfik/dr 
from Eq. 10 are used, and is the time-derivative 
correlation function 0, la 1 a-nd has a natural scaling 
form given by 



Tfc 



d(f)k d<j)^k 
dr dr 
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where L = Lot^^"^, h{x) = C /x is the 2D scaling function 
as X ^ 1, and Lq and C are constants. We can then 
solve for AE in Eq. Q and for the analytic time step 
Ar: 



Ar 



LlAt 
67rSi/3 



dx 



h{x) 



X H- At(ai - l)x^/L^ 
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where x = kL, D = (ai — 1)/Lq and £ 
Solving the integral, we obtain that, 



(11) 
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where ( = LqC^/gi - 1/{12B^/^). The above formula 
is the central result of this manuscript, and implies that 
At < At in general. We now explore how to use this 
result to obtain an accelerated algorithm. 

Scaling of field derivatives in Fourier space: In or- 
der to explore the accuracy of accelerated algorithms 
in Fourier space, it is necessary to know the scaling of 
field derivatives both in the bulk (where k ~ 1/i) and 
near the interface (where k ^ 1 /O • The structure factor 
S{k) = (|0fc|2) = L'^gikL), where g{kL) ~ 1 as fc ~ 1/L 
and g{kL) - (kL)-^ ~ L'^ as A: ~ 1/^0. Therefore we 
obtain 



as k • 
as k • 
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Previous studies d showed that d(j)k/dT = 
{dL/dT)k(j)k as kL ^ 1, so we obtain the form 
for the time-derivative correlation function T(k) — 
{\dcl}k/dt\^) = idL/dTfP{\cj)k\^) = {dL/dTfhi{kL), 
where the scaling function hi{kL) = k^L'^g{kL) ^ 1 as 
k - 1/L, and hi{kL) ~ (fcL)-i ~ L^^ as fc ~ 1/^. 
Therefore we obtain 



dr 



r-2/3 



as k 
as k 



l/L 
1/C 



higher order 

is (|5>fc/ar"|2) 



(14) 
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The generalization of 
derivative correlations 
(dL/dr)2fc2(|9"-i0/c/9r"i|2), where indicates 
that generally the left hand side may not exactly be 
equal to the right hand side. Applying this relation will 
yield (|9"^fc/9r"|2) - {dL/dTf''L^-^''h„{kL), where 
hn{kL) = k'^L'^hn-iikL) ~ {kLf"^'^ - 1 as fc ~ 1/L, 
and hn{kL) ~ (fcL)2"-3 ~ L^^-s as fc ~ 1/^. Therefore 
we have 



Qn 



-ri+l/3 
-2ji/3-1/6 



as k 
as k 



1/L 
1/e 
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The above expression is valid for n > for conserved two 
dimensional scalar order parameter(s). 

Determine the driving scheme for arbitrary accuracy: 
Next, we determine the driving scheme for arbitrary ac- 
curacy in terms of the Fourier space error. Before we 
study the error, we must first distinguish the error in the 
bulk and the error near the interface. Eq. ITJ implies 
Ateff ~ T^/3 as 1/L and Ateff ~ const, as fc ~ 1/^, 
we obtain that the ratio of the single step field update 
with respect to the field At/)^/^^ ~ {At^f fdcj^u / dr) / 4>k is 
of order 0{t-^I^) &skr^l/L and 0(t-2/3) as fc ~ 1/^. 



Therefore the error near the interface is negligible com- 
pared with the error in the bulk, and we will only study 
the error of those modes where k ^ 1/L. 

In Fourier space, we compare the field evolved by an 
unconditionally stable algorithm to the exact dynamics 
evolved by the same amount of energy. Using this crite- 
rion we obtain the Fourier space single step error 

A(j)l = + ^t) ~ (f>k{t + At) 



[Ateff - At) — - 



1 



1 + DP 



(16) 



where Eq. H12|l and d(j)k/dT ~ r"^/^ as fc ~ 1/L are 
used. The values of C and D are finite. Assuming the al- 
gorithmic time step At = At^ , then i = ^/At^/^^^Z^. In 
order to obtain arbitrary accuracy for A0| at arbitrarily 
large t, we require that P — 2/3 since /? > 2/3 will make 
the error uncontrolled (arbitrarily large) at arbitrarily 
large t, and /3 < 2/3 will make the algorithm wastefuUy 
accurate (error is always zero) at arbitrarily large t. A 
is then selected so that a desired accuracy is obtained. 
Thus, At = At^/^ and A0^ - 0{£^) - 0(^^/2). 

For smaU A, Eq. imphes that t t (^1 - CVA^ ■ 
Therefore we can express the algorithmic time step At in 
terms of algorithmic time t: 



At = A(l- CVI - At^/' 



2/3 



(17) 



Writing the driving algorithmic time step in terms of al- 
gorithmic time t instead of the analytic time r has the 
computational advantage of avoiding an intermediate cal- 
culation of T at each update, and thus makes the com- 
putational implementation more straightforward. 

Accuracy in correlations: Lastly, we analytically con- 
firm the numerical results in a previous study |^ that the 
error in structure factor scales as \/A. The Fourier space 
single-step error Eq. I|16|) will at worst accumulate with 
each update. For a small A, evolving to r with time-step 
Ai = v4i^/3 ~ ~ ^j- ^ dr/dn requires a number 

of steps 



= dn ■ 



dr 
A^ 



3t1/3 
A ' 



(18) 



Therefore, at r, we obtain the upper bound on the 
Fourier space multi-step error: 



AC - ^<l>> 



3^3/2^1/3 

A 



lVa. 



(19) 



We can use this to bound the error in the scaled structure 
factor g{kL) — (|0fcp)/L^. As was investigated in our 
previous numerical studies , this quantity is simply the 
magnitude of the difference between the structure factor 
obtained using an unconditionally stable algorithm with 
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one using the exact dynamics at the same energy. As 
k ^ 1/L (in the bulk), we obtain the maximum error: 

Agrna. ~ }\ - VA, (20) 

where (pk ^ L a.s k ^ 1/L is used. Eq. (|20|l is precisely 
the same as the results obtained in our previous, solely 
numerical, study 'Bj. Thus, the error produced in the 
bulk dominates the total error, as it decays much slower 
than the error produced near the interface. This error ac- 
cumulates over time and results the error in the structure 
factor scaling as ^/A. 

In summary, we have analyzed numerical methods for 
solving the Cahn-Hilliard Equation. By explicitly distin- 
guishing the analytic and algorithmic time steps, we have 
developed a relation between them, and have obtained 



an optimal driving scheme At = At^^^ under the require- 
ment of arbitrary accuracy. With this driving scheme, we 
have proved that the upper bound of the multi-step er- 
ror in structure factor scales as \/^, a result obtained by 
numerical methods in a previous study j^] . We note that 
the argument developed herein is founded, ultimately, on 
the physical relationship between the domain size and 
the analytic time (based itself on the scaling hypothe- 
sis). For systems where such relationships exist, or can 
be derived, we expect that this analysis should generalize 
to other systems, such as newly developed Phase Field 
Crystal model 0| . We hope to report this work in a 
subsequent paper. 

MC would like to acknowledge Andrew Rutenberg for 
valuable discussions on this work. 
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